MODULE MINIMIZATION
  USE NRUTIL, ONLY : vabs,assert_eq,nrerror,NR_EPS
  IMPLICIT NONE
  PRIVATE 
  PUBLIC BFGS,FRPR,BFGS_STATA,brent
  
CONTAINS
   
  SUBROUTINE BFGS(theta,gtol,told,func,id,forward)
    IMPLICIT NONE
    REAL(8), INTENT(IN)    :: gtol,told
    INTEGER, INTENT(IN)    :: id        ! For MPI
    REAL(8), INTENT(INOUT) :: theta(:)	! parameter vector over which you which to maximize
    LOGICAL, INTENT(IN) :: forward
    INTERFACE	                        ! This lets the subroutine know that func is not a variable, it is a function
       REAL(8) FUNCTION func(theta)   ! This is not the likelihood itself, it just let the program know
         IMPLICIT NONE		  
         REAL(8), INTENT(IN) :: theta(:)
       END FUNCTION func
    END INTERFACE
    INTEGER, PARAMETER :: ITMAX=4000	! Maximum number of iterations (in case it does not converge)
    REAL(8), PARAMETER :: STPMX=0.1d0,EPS=2.220446049250313D-9 !,TOLD=1.0d-12 ! Tolerance criterions
    INTEGER :: j,its ! indices
    LOGICAL :: check 
    REAL(8) :: den,fac,fad,fae,ftheta,sumdg,sumd ! working matrices
    REAL(8) :: dg(SIZE(theta)),g(SIZE(theta)),hdg(SIZE(theta)),thetanew(SIZE(theta)),d(SIZE(theta))
    REAL(8) :: B(SIZE(theta),SIZE(theta)),fret,stpmax,tim
    INTEGER :: tin,tout,hz,n
    CHARACTER(len=500) :: file_open
    n=SIZE(theta)
    ftheta=func(theta)        ! initial value of the function
    IF (id==0) WRITE(*,'(A,F16.8)') 'First value of the likelihood ',ftheta
    CALL FLUSH(6)
    IF (forward) THEN
       g=gradient1(theta,ftheta,func)	! initial value of the gradient, if analytical gradient known here
    ELSE
       g=gradient2(theta,func)	! initial value of the gradient, if analytical gradient known here
    END IF
    B=0.0d0				! initialize the matrix to a negative definite matrix
    DO j = 1, n
       B(j,j)=1.0d0
       IF (id==0) WRITE(*,*) j,g(j)
       CALL FLUSH(6)
    END DO
    d=-MATMUL(B,g)	! initial line direction
    stpmax=STPMX*MAX(vabs(theta),DBLE(n))
    DO its=1,ITMAX	! Main loop over iterations
       CALL SYSTEM_CLOCK(count_rate=hz) 
       CALL SYSTEM_CLOCK(count=tin) 
       CALL lnsrch(theta,ftheta,g,d,thetanew,fret,stpmax,check,func)
       ftheta=fret	! update the value of the fuction
       d=thetanew-theta	! update the direction of movement
       theta=thetanew	! update the point
       IF (MAXVAL(ABS(d)/MAX(ABS(theta),1.0D0)) < TOLD) THEN ! Check convergence in Delta-d
          IF (id==0) THEN
             WRITE(*,*) 'TOLD'
             OPEN(6424,file="point.out")
             OPEN(6425,file="hess.out")
             DO j = 1, n
                WRITE(6424,'(F32.16)') theta(j)
                WRITE(6425,'(2000F32.16)') B(j,j)
             END DO
             CLOSE(6424)	
             CLOSE(6425)	
          END IF
          RETURN     
       END IF
       dg=g				! save the old gradient
       IF (forward) THEN
          g=gradient1(theta,ftheta,func)	! and get the new gradient
       ELSE
          g=gradient2(theta,func)	! and get the new gradient
       END IF
       den=MAX(ftheta,1.0d0)	
       IF (MAXVAL(ABS(g)*MAX(ABS(theta),1.0d0)/den) < gtol) THEN ! Check for convergence on zero gradient
          IF (id==0) THEN
             WRITE(*,*) 'GTOL'
             OPEN(6424,file="point.out")
             OPEN(6425,file="hess.out")
             DO j = 1, n
                WRITE(6424,'(F32.16)') theta(j)
                WRITE(6425,'(2000F32.16)') B(j,j)
             END DO
             CLOSE(6424)	
             CLOSE(6425)	
          END IF
          RETURN
       END IF
       dg=g-dg			! compute difference in gradients
       hdg=MATMUL(B,dg)	        ! and difference times current matrix
       fac=DOT_PRODUCT(dg,d)	! Dot products for denominators
       fae=DOT_PRODUCT(dg,hdg)
       sumdg=DOT_PRODUCT(dg,dg)
       sumd=DOT_PRODUCT(d,d)
       IF (fac > SQRT(EPS*sumdg*sumd)) THEN	! Skip update if fac not positive enough
          fac=1.0d0/fac
          fad=1.0d0/fae
          dg=fac*d-fad*hdg			! Vector that makes BFGS different from DFP
          ! BFGS updating formula
          B=B+fac*P_Outer_Product_Min(d,d)
          B=B-fad*P_Outer_Product_Min(hdg,hdg)
          B=B+fae*P_Outer_Product_Min(dg,dg)
       END IF
       d=-MATMUL(B,g)	! Next direction to go
       IF (id==0) THEN
          OPEN(6424,file="point.out")
          DO j = 1, n
             WRITE(6424,'(F32.16)') theta(j)
          END DO
          CLOSE(6424)	
       END IF
       CALL SYSTEM_CLOCK(count=tout) 
       tout=tout-tin
       tim=DBLE(tout)/DBLE(hz)
       IF (id==0) THEN
          write(*,'(A,I10,A,F32.16,A,F16.6)') 'iteration:',its,' value:',fret,' time:',tim
          CALL FLUSH(6)
       END IF
    END DO		! Go back for another iteration
!    CALL NRERROR('BFGS: too many iterations')
    IF (id==0) THEN
       WRITE(*,*) 'GTOL'
       OPEN(6424,file="point.out")
       OPEN(6425,file="hess.out")
       DO j = 1, n
          WRITE(6424,'(F32.16)') theta(j)
          WRITE(6425,'(2000F32.16)') B(j,j)
       END DO
       CLOSE(6424)	
       CLOSE(6425)	
    END IF
    WRITE(*,*) 'BFGS:Too many iterations'
  END SUBROUTINE BFGS
  
  SUBROUTINE BFGS_STATA(theta,gtol,told,func,id,forward)
    IMPLICIT NONE
    REAL(8), INTENT(IN)    :: gtol,told
    INTEGER, INTENT(IN)    :: id        ! For MPI
    REAL(8), INTENT(INOUT) :: theta(:)	! parameter vector over which you which to maximize
    LOGICAL, INTENT(IN) :: forward
    INTERFACE	                        ! This lets the subroutine know that func is not a variable, it is a function
       REAL(8) FUNCTION func(theta)   ! This is not the likelihood itself, it just let the program know
         IMPLICIT NONE		  
         REAL(8), INTENT(IN) :: theta(:)
       END FUNCTION func
    END INTERFACE
    INTEGER, PARAMETER :: ITMAX=100	! Maximum number of iterations (in case it does not converge)
    REAL(8), PARAMETER :: STPMX=0.01d0,EPS=2.220446049250313D-9 !,TOLD=1.0d-12 ! Tolerance criterions
    INTEGER :: j,its ! indices
    LOGICAL :: check 
    REAL(8) :: den,fac,fad,fae,ftheta,sumdg,sumd ! working matrices
    REAL(8) :: dg(SIZE(theta)),g(SIZE(theta)),hdg(SIZE(theta)),thetanew(SIZE(theta)),d(SIZE(theta))
    REAL(8) :: B(SIZE(theta),SIZE(theta)),fret,stpmax,tim,SD(SIZE(theta))
    INTEGER :: tin,tout,hz,n
    CHARACTER(len=500) :: file_open
    n=SIZE(theta)
    ftheta=func(theta)        ! initial value of the function
    IF (id==0) WRITE(*,'(A,F16.8)') 'First value of the likelihood ',ftheta
    CALL FLUSH(6)
    IF (forward) THEN
       g=gradient1(theta,ftheta,func)	! initial value of the gradient, if analytical gradient known here
    ELSE
       g=gradient2(theta,func)	! initial value of the gradient, if analytical gradient known here
    END IF
    B=0.0d0				! initialize the matrix to a negative definite matrix
    DO j = 1, n
       B(j,j)=1.0d0
       IF (id==0) WRITE(*,*) j,g(j)
       CALL FLUSH(6)
    END DO
    d=-MATMUL(B,g)	! initial line direction
    stpmax=STPMX*MAX(vabs(theta),DBLE(n))
    DO its=1,ITMAX	! Main loop over iterations
       CALL SYSTEM_CLOCK(count_rate=hz) 
       CALL SYSTEM_CLOCK(count=tin) 
       CALL lnsrch(theta,ftheta,g,d,thetanew,fret,stpmax,check,func)
       ftheta=fret	! update the value of the fuction
       d=thetanew-theta	! update the direction of movement
       theta=thetanew	! update the point
       IF (MAXVAL(ABS(d)/MAX(ABS(theta),1.0D0)) < TOLD) THEN ! Check convergence in Delta-d
          IF (id==0) THEN
             WRITE(*,*) 'TOLD'
             OPEN(6424,file="point.out")
             DO j = 1, n
                WRITE(6424,'(F32.16)') theta(j)
             END DO
             CLOSE(6424)	
          END IF
          RETURN     
       END IF
       dg=g				! save the old gradient
       IF (forward) THEN
          g=gradient1(theta,ftheta,func)	! and get the new gradient
       ELSE
          g=gradient2(theta,func)	! and get the new gradient
       END IF
       den=MAX(ftheta,1.0d0)	
       IF (MAXVAL(ABS(g)*MAX(ABS(theta),1.0d0)/den) < gtol) THEN ! Check for convergence on zero gradient
          IF (id==0) THEN
             WRITE(*,*) 'GTOL'
             OPEN(6424,file="point.out")
             DO j = 1, n
                WRITE(6424,'(F32.16)') theta(j)
             END DO
             CLOSE(6424)	
          END IF
          RETURN
       END IF
       dg=g-dg			! compute difference in gradients
       hdg=MATMUL(B,dg)	        ! and difference times current matrix
       fac=DOT_PRODUCT(dg,d)	! Dot products for denominators
       fae=DOT_PRODUCT(dg,hdg)
       sumdg=DOT_PRODUCT(dg,dg)
       sumd=DOT_PRODUCT(d,d)
       IF (ABS(fac*fae) > EPS*sumdg*sumd) THEN	! Skip update if fac not positive enough
          fac=1.0d0/fac
          fad=1.0d0/fae
          SD=d*fac - hdg*fad
          B = B - fac*P_Outer_Product_Min(d,d)
          B = B + fae*P_Outer_Product_Min(SD,SD)
          B = B - fad*P_Outer_Product_Min(hdg,hdg)
          B = 0.5d0*(B + TRANSPOSE(B))
       END IF
       d=-MATMUL(B,g)	! Next direction to go
       IF (id==0) THEN
          OPEN(6424,file="point.out")
          DO j = 1, n
             WRITE(6424,'(F32.16)') theta(j)
          END DO
          CLOSE(6424)	
       END IF
       CALL SYSTEM_CLOCK(count=tout) 
       tout=tout-tin
       tim=DBLE(tout)/DBLE(hz)
       IF (id==0) THEN
          write(*,'(A,I10,A,F32.16,A,F16.6)') 'iteration:',its,' value:',fret,' time:',tim
          CALL FLUSH(6)
       END IF
    END DO		! Go back for another iteration
!    CALL NRERROR('BFGS: too many iterations')
    WRITE(*,*) 'BFGS_STATA:Too many iterations'
  END SUBROUTINE BFGS_STATA

  SUBROUTINE FRPR(p,ftol,iter,fret,func,id,forward)
    USE nrutil, ONLY : nrerror
    IMPLICIT NONE
    INTERFACE
       FUNCTION func(p)
         IMPLICIT NONE
         REAL(8), DIMENSION(:), INTENT(IN) :: p
         REAL(8) :: func
       END FUNCTION func
    END INTERFACE
    LOGICAL, INTENT(IN) :: forward
    INTEGER, INTENT(IN) :: id
    INTEGER, INTENT(OUT) :: iter
    REAL(8), INTENT(IN) :: ftol
    REAL(8), INTENT(OUT) :: fret
    REAL(8), DIMENSION(:), INTENT(INOUT) :: p
    INTEGER, PARAMETER :: ITMAX=200
    REAL(8), PARAMETER :: EPS=1.0d-10
    INTEGER :: its,tin,tout,hz,j
    REAL(8) :: dgg,fp,gam,gg,tim
    REAL(8), DIMENSION(size(p)) :: g,h,xi
    fp=func(p)
    IF (id==0) WRITE(*,'(A,F16.8)') 'First value of the likelihood ',fp
    CALL FLUSH(6)
    IF (forward) THEN
       xi=gradient1(p,fp,func)	! initial value of the gradient, if analytical gradient known here
    ELSE
       xi=gradient2(p,func)
    END IF
    DO j = 1,SIZE(p)
       IF (id==0) WRITE(*,*) j,xi(j)
       CALL FLUSH(6)
    END DO
    g=-xi
    h=g
    xi=h
    DO its=1,ITMAX
       CALL SYSTEM_CLOCK(count_rate=hz) 
       CALL SYSTEM_CLOCK(count=tin) 
       iter=its
       CALL LINMIN(p,xi,fret,func)
       IF (id==0) THEN
          WRITE(*,*) 'Linmin'
          CALL FLUSH(6)
          OPEN(6424,file="point.out")
          DO j = 1, SIZE(p)
             WRITE(6424,'(F32.16)') p(j)
          END DO
          CLOSE(6424)	
       END IF
       IF (2.0d0*abs(fret-fp) <= ftol*(ABS(fret)+ABS(fp)+EPS)) RETURN
       fp=fret
       IF (forward) THEN
          xi=gradient1(p,fp,func)	! initial value of the gradient, if analytical gradient known here
       ELSE
          xi=gradient2(p,func)
       END IF
       gg=DOT_PRODUCT(g,g)
       !    dgg=dot_product(xi,xi)
       dgg=DOT_PRODUCT(xi+g,xi)
       IF (gg == 0.0d0) RETURN
       gam=dgg/gg
       g=-xi
       h=g+gam*h
       xi=h
       CALL SYSTEM_CLOCK(count=tout) 
       tout=tout-tin
       tim=DBLE(tout)/DBLE(hz)
       IF (id==0) THEN
          write(*,'(A,I10,A,F32.16,A,F16.6)') 'iteration:',its,' value:',fret,' time:',tim
          CALL FLUSH(6)
       END IF
    END DO
    CALL nrerror('frpr: maximum iterations exceeded')
  END SUBROUTINE FRPR
  
  SUBROUTINE linmin(p,xi,fret,func)
    USE nrutil, ONLY : assert_eq
    USE FRPR_PRIVATE_F1DIM
    IMPLICIT NONE
    INTERFACE
       FUNCTION func(p)
         IMPLICIT NONE
         REAL(8), DIMENSION(:), INTENT(IN) :: p
         REAL(8) :: func
       END FUNCTION func
    END INTERFACE
    REAL(8), INTENT(OUT) :: fret
    REAL(8), DIMENSION(:), TARGET, INTENT(INOUT) :: p,xi
    REAL(8), PARAMETER :: TOL=1.0d-4
    REAL(8) :: ax,bx,fa,fb,fx,xmin,xx
    ncom=assert_eq(size(p),size(xi),'linmin')
    pcom=>p
    xicom=>xi
    ax=0.0d0
    xx=1.0d0
    CALL mnbrak(ax,xx,bx,fa,fx,fb,f1dim,func)
    WRITE(*,*) 'mnbrak'
    fret=brent(ax,xx,bx,f1dim,TOL,xmin,func)
    xi=xmin*xi
    p=p+xi
  END SUBROUTINE linmin
  
  FUNCTION brent(ax,bx,cx,feval,tol,xmin,func)
    USE nrutil, ONLY : nrerror
    IMPLICIT NONE
    REAL(8), INTENT(IN) :: ax,bx,cx,tol
    REAL(8), INTENT(OUT) :: xmin
    REAL(8) :: brent
    INTERFACE
       FUNCTION func(p)
         IMPLICIT NONE
         REAL(8), DIMENSION(:), INTENT(IN) :: p
         REAL(8) :: func
       END FUNCTION func
       FUNCTION feval(x,func)
         IMPLICIT NONE
         INTERFACE
            FUNCTION func(p)
              IMPLICIT NONE
              REAL(8), DIMENSION(:), INTENT(IN) :: p
              REAL(8) :: func
            END FUNCTION func
         END INTERFACE
         REAL(8), INTENT(IN) :: x
         REAL(8) :: feval
       END FUNCTION feval
    END INTERFACE
    INTEGER, PARAMETER :: ITMAX=100
    REAL(8), PARAMETER :: CGOLD=0.3819660d0,ZEPS=1.0d-3*epsilon(ax)
    INTEGER :: iter
    REAL(8) :: a,b,d,e,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm
    a=MIN(ax,cx)
    b=MAX(ax,cx)
    v=bx
    w=v
    x=v
    e=0.0d0
    fx=feval(x,func)
    fv=fx
    fw=fx
    DO iter=1,ITMAX
       xm=0.5d0*(a+b)
       tol1=tol*ABS(x)+ZEPS
       tol2=2.0d0*tol1
       IF (ABS(x-xm) <= (tol2-0.5d0*(b-a))) THEN
          xmin=x
          brent=fx
          RETURN
       END IF
       IF (ABS(e) > tol1) THEN
          r=(x-w)*(fx-fv)
          q=(x-v)*(fx-fw)
          p=(x-v)*q-(x-w)*r
          q=2.0d0*(q-r)
          IF (q > 0.0d0) p=-p
          q=ABS(q)
          etemp=e
          e=d
          IF (ABS(p) >= ABS(0.5d0*q*etemp) .OR. &
               p <= q*(a-x) .OR. p >= q*(b-x)) THEN
             e=MERGE(a-x,b-x, x >= xm )
             d=CGOLD*e
          ELSE
             d=p/q
             u=x+d
             IF (u-a < tol2 .OR. b-u < tol2) d=SIGN(tol1,xm-x)
          END IF
       ELSE
          e=MERGE(a-x,b-x, x >= xm )
          d=CGOLD*e
       END IF
       u=MERGE(x+d,x+SIGN(tol1,d), ABS(d) >= tol1 )
       fu=feval(u,func)
       IF (fu <= fx) THEN
          IF (u >= x) THEN
             a=x
          ELSE
             b=x
          END IF
          CALL shft(v,w,x,u)
          CALL shft(fv,fw,fx,fu)
       ELSE
          IF (u < x) THEN
             a=u
          ELSE
             b=u
          END IF
          IF (fu <= fw .OR. w == x) THEN
             v=w
             fv=fw
             w=u
             fw=fu
          ELSE IF (fu <= fv .OR. v == x .OR. v == w) THEN
             v=u
             fv=fu
          END IF
       END IF
    END DO
    CALL nrerror('brent: exceed maximum iterations')
  CONTAINS
    
    SUBROUTINE shft(a,b,c,d)
      REAL(8), INTENT(OUT) :: a
      REAL(8), INTENT(INOUT) :: b,c
      REAL(8), INTENT(IN) :: d
      a=b
      b=c
      c=d
    END SUBROUTINE shft
  END FUNCTION brent
  
  SUBROUTINE mnbrak(ax,bx,cx,fa,fb,fc,feval,func)
    USE nrutil, ONLY : swap
    IMPLICIT NONE
    REAL(8), INTENT(INOUT) :: ax,bx
    REAL(8), INTENT(OUT) :: cx,fa,fb,fc
    INTERFACE
       FUNCTION func(p)
         IMPLICIT NONE
         REAL(8), DIMENSION(:), INTENT(IN) :: p
         REAL(8) :: func
       END FUNCTION func
       FUNCTION feval(x,func)
         IMPLICIT NONE
         INTERFACE
            FUNCTION func(p)
              IMPLICIT NONE
              REAL(8), DIMENSION(:), INTENT(IN) :: p
              REAL(8) :: func
            END FUNCTION func
         END INTERFACE
         REAL(8), INTENT(IN) :: x
         REAL(8) :: feval
       END FUNCTION feval
    END INTERFACE
    REAL(8), PARAMETER :: GOLD=1.618034d0,GLIMIT=100.0d0,TINY=1.0d-20
    REAL(8) :: fu,q,r,u,ulim
    fa=feval(ax,func)
    fb=feval(bx,func)
    IF (fb > fa) THEN
       CALL swap(ax,bx)
       CALL swap(fa,fb)
    END IF
    cx=bx+GOLD*(bx-ax)
    fc=feval(cx,func)
    DO
       IF (fb < fc) RETURN
       r=(bx-ax)*(fb-fc)
       q=(bx-cx)*(fb-fa)
       u=bx-((bx-cx)*q-(bx-ax)*r)/(2.0d0*sign(max(abs(q-r),TINY),q-r))
       ulim=bx+GLIMIT*(cx-bx)
       if ((bx-u)*(u-cx) > 0.0d0) then
          fu=feval(u,func)
          if (fu < fc) then
             ax=bx
             fa=fb
             bx=u
             fb=fu
             RETURN
          else if (fu > fb) then
             cx=u
             fc=fu
             RETURN
          end if
          u=cx+GOLD*(cx-bx)
          fu=feval(u,func)
       else if ((cx-u)*(u-ulim) > 0.0d0) then
          fu=feval(u,func)
          if (fu < fc) then
             bx=cx
             cx=u
             u=cx+GOLD*(cx-bx)
             call shft(fb,fc,fu,feval(u,func))
          end if
       else if ((u-ulim)*(ulim-cx) >= 0.0d0) then
          u=ulim
          fu=feval(u,func)
       else
          u=cx+GOLD*(cx-bx)
          fu=feval(u,func)
       end if
       call shft(ax,bx,cx,u)
       call shft(fa,fb,fc,fu)
    end do
  CONTAINS
    
    SUBROUTINE shft(a,b,c,d)
      REAL(8), INTENT(OUT) :: a
      REAL(8), INTENT(INOUT) :: b,c
      REAL(8), INTENT(IN) :: d
      a=b
      b=c
      c=d
    END SUBROUTINE shft
  END SUBROUTINE mnbrak

  FUNCTION GRADIENT1(x0,f0,func)
    ! Routine to get the numerical gradient of a function using forward differences
    IMPLICIT NONE
    REAL(8), INTENT(IN) :: x0(:),f0
    INTERFACE	                       ! This lets the subroutine know that func is not a variable, it is a function
       REAL(8) FUNCTION func(theta)  ! This is not the likelihood itself, it just let the program know
         IMPLICIT NONE		       ! the arguments and type of arguments it accepts
         REAL(8), INTENT(IN) :: theta(:)
       END FUNCTION func
    END INTERFACE
    REAL(8), PARAMETER :: EPS=1.49011611938477D-08 !EPS=SQRT(NR_EPS)
    REAL(8) :: GRADIENT1(SIZE(x0)),grdd(SIZE(x0)),dh(SIZE(x0)),xdh(SIZE(x0))
    REAL(8) :: arg(SIZE(x0),SIZE(x0))
    INTEGER :: i,k
    k=SIZE(x0)
    grdd = 0.0D0
    ! Computation of stepsize (dh) for gradient
    dh = EPS*ABS(x0)
    WHERE (dh<1.0d-13) dh=EPS
    xdh = x0+dh
    dh=xdh-x0
    arg = SPREAD(x0,DIM=2,NCOPIES=k)
    DO i = 1, k
       arg(i,i)=xdh(i)    
       grdd(i)=func(arg(:,i))
    END DO
    GRADIENT1 = (grdd-f0)/dh
  END FUNCTION GRADIENT1
  
  FUNCTION GRADIENT2(x0,func)
    ! Routine to get the numerical gradient of a function using central differences
    IMPLICIT NONE
    REAL(8), INTENT(IN) :: x0(:)
    INTERFACE	                       ! This lets the subroutine know that func is not a variable, it is a function
       REAL(8) FUNCTION func(theta)  ! This is not the likelihood itself, it just let the program know
         IMPLICIT NONE		       ! the arguments and type of arguments it accepts
         REAL(8), INTENT(IN) :: theta(:)
       END FUNCTION func
    END INTERFACE
    REAL(8), PARAMETER :: EPS=6.05545445239334D-06 !EPS=NR_EPS**(1.0d0/3.0d0)
    REAL(8) :: GRADIENT2(SIZE(x0)),grdd(SIZE(x0)),dh(SIZE(x0)),xdh1(SIZE(x0)),xdh0(SIZE(x0)),dhh(SIZE(x0))
    REAL(8) :: arg(SIZE(x0),SIZE(x0))
    INTEGER :: i,k
    k=SIZE(x0)
    grdd = 0.0D0
    ! Computation of stepsize (dh) for gradient
    dh = EPS*ABS(x0)
    WHERE (dh==0.0D0) dh=EPS
    xdh1 = x0+dh
    arg = SPREAD(x0,DIM=2,NCOPIES=k)
    DO i = 1, k
       arg(i,i)=xdh1(i)    
       grdd(i)=func(arg(:,i))
    END DO
    xdh0 = x0-dh
    dhh = xdh1-xdh0
    arg = SPREAD(x0,DIM=2,NCOPIES=k)
    DO i = 1, k
       arg(i,i)=xdh0(i)    
       grdd(i)=grdd(i)-func(arg(:,i))
    END DO
    GRADIENT2 = grdd/dhh
  END FUNCTION GRADIENT2

  FUNCTION GRADIENT3(x0,func)
    ! Routine to get the numerical gradient of a function using a 5 point stencil
    IMPLICIT NONE
    REAL(8), INTENT(IN) :: x0(:)
    INTERFACE	                       ! This lets the subroutine know that func is not a variable, it is a function
       REAL(8) FUNCTION func(theta)  ! This is not the likelihood itself, it just let the program know
         IMPLICIT NONE		       ! the arguments and type of arguments it accepts
         REAL(8), INTENT(IN) :: theta(:)
       END FUNCTION func
    END INTERFACE
    REAL(8), PARAMETER :: EPS=6.05545445239334D-06 !EPS=NR_EPS**(1.0d0/3.0d0)
    REAL(8) :: GRADIENT3(SIZE(x0)),grdd(SIZE(x0)),dh(SIZE(x0)),xdh(SIZE(x0))
    REAL(8) :: arg(SIZE(x0),SIZE(x0))
    INTEGER :: i,k
    k=SIZE(x0)
    grdd = 0.0D0
    ! Computation of stepsize (dh) for gradient
    dh = EPS*ABS(x0)
    WHERE (dh==0.0D0) dh=EPS
    xdh = x0+dh
    dh = xdh-x0
    arg = SPREAD(x0,DIM=2,NCOPIES=k)
    DO i = 1, k
       arg(i,i)=xdh(i)    
       grdd(i)=func(arg(:,i))
    END DO
    xdh = x0-dh
    dh = x0-xdh
    arg = SPREAD(x0,DIM=2,NCOPIES=k)
    DO i = 1, k
       arg(i,i)=xdh(i)    
       grdd(i)=grdd(i)-func(arg(:,i))
    END DO
    GRADIENT3 = grdd/(2.0d0*dh)
  END FUNCTION GRADIENT3

  SUBROUTINE lnsrch(xold,fold,g,p,x,f,stpmax,check,func)
    ! Given an N-dimensional point xold, the value of the function and gradient there, fold
    ! and g, and a direction p, finds a new point x along the direction p from xold where the
    ! function func has decreased sufficiently. xold, g, p, and x are all arrays of length N.
    ! The new function value is returned in f. stpmax is an input quantity that limits the length
    ! of the steps so that you do not try to evaluate the function in regions where it is undefined
    ! or subject to overflow. p is usually the Newton direction. The output quantity check is
    ! false on a normal exit. It is true when x is too close to xold. In a minimization algorithm,
    ! this usually signals convergence and can be ignored. However, in a zero-finding algorithm
    ! the calling program should check whether the convergence is spurious.
    ! Parameters: ALF ensures su.cient decrease in function value; TOLX is the convergence
    ! criterion on .x.
    IMPLICIT NONE
    REAL(8), DIMENSION(:), INTENT(IN) :: xold,g
    REAL(8), DIMENSION(:), INTENT(INOUT) :: p(SIZE(xold))
    REAL(8), INTENT(IN) :: fold,stpmax
    REAL(8), INTENT(OUT) :: x(SIZE(xold))
    REAL(8), INTENT(OUT) :: f
    LOGICAL, INTENT(OUT) :: check
    INTERFACE	! This lets the subroutine know that func is not a variable, it is a function
       REAL(8) FUNCTION func(theta)  ! This is not the likelihood itself, it just let the program know
         IMPLICIT NONE		
         REAL(8), INTENT(IN) :: theta(:)
       END FUNCTION func
    END INTERFACE
    REAL(8), PARAMETER :: ALF=1.0d-4,TOLX=1.0d-9
    INTEGER :: ndum
    REAL(8) :: a,alam,alam2,alamin,b,disc,f2,pabs,rhs1,rhs2,slope,tmplam,zero
    zero=0.0D0
    ndum=assert_eq(SIZE(g),SIZE(xold),'Dimensions do not agree: LNSRCH')
    check=.false.
    pabs=vabs(p(:))
    IF (pabs > stpmax) p(:)=p(:)*stpmax/pabs ! Scale if attempted step is too big.
    slope=DOT_PRODUCT(g,p)
    IF (slope >= zero) CALL NRERROR('roundoff problem in lnsrch')
    alamin=TOLX/MAXVAL(DABS(p(:))/MAX(DABS(xold(:)),1.0D0)) !Compute .min.
    alam=1.0D0 ! Always try full Newton step first.
    DO ! Start of iteration loop.
       x(:)=xold(:)+alam*p(:)
       F=FUNC(x)
       IF (alam < alamin) THEN                 ! Convergence on .x. For zero finding,
          !the calling program should verify the convergence.
          x(:)=xold(:)
          check=.true.
          RETURN
       ELSE IF (f <= fold+ALF*alam*slope) THEN ! Sufficient function decrease.
          RETURN
       ELSE !Backtrack.
          IF (alam == 1.0D0) THEN !First time.
             tmplam=-slope/(2.0D0*(f-fold-slope))
          ELSE !Subsequent backtracks.
             rhs1=f-fold-alam*slope
             rhs2=f2-fold-alam2*slope
             a=(rhs1/alam**2-rhs2/alam2**2)/(alam-alam2)
             b=(-alam2*rhs1/alam**2+alam*rhs2/alam2**2)/&
                  (alam-alam2)
             IF (a == 0.0D0) THEN
                tmplam=-slope/(2.0D0*b)
             ELSE
                disc=b*b-3.0*a*slope
                IF (disc < 0.0D0) THEN
                   tmplam=0.5D0*alam
                ELSE IF (b <= 0.0D0) THEN
                   tmplam=(-b+SQRT(disc))/(3.0D0*a)
                ELSE
                   tmplam=-slope/(b+SQRT(disc))
                END IF
             END IF
             IF (tmplam > 0.5D0*alam) tmplam=0.5D0*alam
          END IF
       END IF
       alam2=alam
       f2=f
       alam=MAX(tmplam,0.1D0*alam) 
    END DO !Try again.
  END SUBROUTINE lnsrch
  
  FUNCTION P_Outer_Product_Min(a,b)
    IMPLICIT NONE
    REAL(8), DIMENSION(:), INTENT(IN) :: a,b
    REAL(8), DIMENSION(size(a),size(b)) :: P_Outer_Product_Min
    P_Outer_Product_Min = spread(a,dim=2,ncopies=size(b)) * &
         spread(b,dim=1,ncopies=size(a))
  END FUNCTION P_Outer_Product_Min
  
END MODULE MINIMIZATION
